#include "2BWTSearch.h"
#include "2BWTSearch_cpu.h"
#include "cuda_header.h"
#include <algorithm>
#include "coproc.h"
#include "../CompressAlign/cpp/compressTitle_cpu.h"
#include "../CompressAlign/cpp/compressOutput_cpu.h"
#include "../CompressAlign/cpp/compressOutput_gpu.h"
#include "../saIndex/alignerOutput.h"
#include "result.h"
#include <mian_util.h>
#include <AlignResult.h>
#include "alignerOutput.h"
#include <algorithm>

//#include "Index.h"
//#include "PausedQueryBuffer.h"


using namespace std;

static double inputSec = 0.0;
static double alignSec = 0.0;
static double comprSec = 0.0;
static double outputSec = 0.0;

/*
char getComplement(char base) {

	switch(base) {
		case 'a':
			return 't';
		case 'A':
			return 'T';
		case 't':
			return 'a';
		case 'T':
			return 'A';
		case 'c':
			return 'g';
		case 'C':
			return 'G';
		case 'g':
			return 'c';
		case 'G':
			return 'C';
		default: {
				 printf("character : %c \n", base);
				 ERROR_EXIT("unknown base");
			 }
	}
}
void reverse_complement(char* outRead, const char* inRead, const int readLength) {
	for(int i = 0; i < readLength; i++) {
		outRead[i] = getComplement(inRead[i]);
	}
	reverse(outRead, outRead + readLength);
}

*/

void copyFileName(char * name, char * path, char * suffix)
{
	strcat(name, path);
	strcat(name,suffix);
}

char decodeInt(int i)
{
	switch(i)
	{
		case(0):
			return 'A';
			break;
		case(1):
			return'C';
			break;
		case(2):
			return'G';
			break;
		case(3):
			return'T';
			break;
	}
}

char EncodeInt(char a)
{
	switch(a)
	{
		case'A':
			return 0;
		case'C':
			return 1;
		case'G':
			return 2;
		case'T':
			return 3;
	}
	return 0;
}


uint64 reportNumMatch(const long* numMatch, const int numQuery) {
	int numDiscard = 0;
	int numAlign = 0;
	int numQueryAligned = 0;

	for(int i = 0; i < numQuery; i++) {
		if(numMatch[i] == -1) {
			numDiscard++;
		} else {
			numAlign += numMatch[i];
			if(numMatch[i] > 0) {
				numQueryAligned ++;
			}
		}
	}

	printf("reported by reportNumMatch: numDiscard: %d\n", numDiscard);
	printf("reported by reportNumMatch: numQueryAligned: %d\n", numQueryAligned);
	printf("reported by reportNumMatch: numAlign: %d\n", numAlign);

	return numAlign;
}

void LoadIndex(BWT * bwt, BWT *revbwt, char* path)
{
	int i;
	FILE *bwtCodeFile,*saValueFile, *revbwtCodeFile, *occValueFile,*revoccValueFile, *gpuOccValueFile, *gpurevOccValueFile= NULL;
	unsigned int tmp;
	unsigned int bwtCodeLengthInFile;
	const int MaxFileNameLength = 100;
	//	char path[] = "/home/lumian/data/human_ref/ref_good.fa.index";
	//	char path[] = "/home/lumian/data/human_ref/2bwt-index/hg19.fa.index";
	//char path[] = "data/ref.fa.index";
	char bwtCodeFileNameSuffix[] = ".bwt";
	//	char saValueFileNameSuffix[] = "/ref.fa.index.sa";
	char revbwtCodeFileNameSuffix[] = ".rev.bwt";
	char lookupTableFileNameSuffix[] = ".lkt";
	char revLookupTableFileNameSuffix[] = ".rev.lkt";
	char occValueFileNameSuffix[] = ".fmv";
	char revoccValueFileNameSuffix[] = ".rev.fmv";
	char gpuoccValueFileNameSuffix[] = ".fmv.gpu";
	char gpurevoccValueFileNameSuffix[] = ".rev.fmv.gpu";
	//	char gpuoccValueFileNameSuffix[] = ".fmv.gpu";
	//	char gpurevoccValueFileNameSuffix[] = ".rev.fmv.gpu";


	char bwtCodeFileName[MaxFileNameLength]="";
	char saValueFileName[MaxFileNameLength]="";
	char revbwtCodeFileName[MaxFileNameLength]="";
	char lookupTableFileName[MaxFileNameLength]="";
	char revLookupTableFileName[MaxFileNameLength]="";
	char occValueFileName[MaxFileNameLength]="";
	char revoccValueFileName[MaxFileNameLength]="";

	//char gpuOccValueFileName[] = "/home/lumian/data/human_ref/hg19_OccValueFileGPU";

	char gpuOccValueFileName[MaxFileNameLength] = "";
	char gpurevOccValueFileName[MaxFileNameLength] = "";

	//char gpuOccValueFileName[] = "OccValueFileGPU";
	//char gpurevOccValueFileName[] = "RevOccValueFileGPU";

	copyFileName(bwtCodeFileName, path, bwtCodeFileNameSuffix);
	//	copyFileName(saValueFileName, path, saValueFileNameSuffix);
	copyFileName(revbwtCodeFileName, path, revbwtCodeFileNameSuffix);
	copyFileName(lookupTableFileName, path, lookupTableFileNameSuffix);
	copyFileName(revLookupTableFileName, path, revLookupTableFileNameSuffix);
	copyFileName(occValueFileName, path, occValueFileNameSuffix);
	copyFileName(revoccValueFileName, path, revoccValueFileNameSuffix);
	copyFileName(gpuOccValueFileName, path, gpuoccValueFileNameSuffix);
	copyFileName(gpurevOccValueFileName, path, gpurevoccValueFileNameSuffix);

	unsigned int *tempBwtCode;
	unsigned int * tempPacValue;

	bwtCodeFile = (FILE*)fopen(bwtCodeFileName, "rb");
	if (bwtCodeFile == NULL) {
		fprintf(stderr, "BWTLoad() : cannot open bwtCodeFile%s!\n", bwtCodeFileName);
		exit(1);
	}

	revbwtCodeFile = fopen(revbwtCodeFileName, "rb");
	if (revbwtCodeFile == NULL) {
		fprintf(stderr, "BWTLoad() : cannot open revbwtCodeFile!\n");
		exit(1);
	}
#if LAYER_INDEX_2
	occValueFile = (FILE*)fopen(occValueFileName, "rb");
	if (occValueFile == NULL) {
		fprintf(stderr, "BWTLoad() : cannot open occFile!\n");
		exit(1);
	}

	revoccValueFile = fopen(revoccValueFileName, "rb");
	if (revoccValueFile == NULL) {
		fprintf(stderr, "BWTLoad() : cannot open revoccFile!\n");
		exit(1);
	}
#endif

#if LAYER_INDEX_1
	gpuOccValueFile = (FILE *)fopen(gpuOccValueFileName, "rb");
	gpurevOccValueFile = (FILE *)fopen(gpurevOccValueFileName,"rb");
#endif
	/*
	   if (saValueFileName != NULL && saValueFileName[0] != '\0' && saValueFileName[0] != '-') {
	   saValueFile = (FILE*)fopen(saValueFileName, "rb");
	   if (saValueFile == NULL) {
	   fprintf(stderr, "BWTLoad() : cannot open saValueFile!\n");
	   fprintf(stderr, "saValueFileName:%s\n", saValueFileName);
	   exit(1);
	   }
	   }

	 */
	printf("reference file: %s\n",path);
	bwt->countValue = (unsigned int *)malloc((ALPHABET_SIZE + 1)*sizeof(int));
	fread(&bwt->inverseSa0, sizeof(unsigned int), 1, bwtCodeFile);
	fread(&bwt->countValue[1], sizeof(unsigned int), ALPHABET_SIZE, bwtCodeFile);
	Initialize(bwt, bwt->countValue[ALPHABET_SIZE]);
	bwt->countValue[0] = 0;

	bwt->bwtSizeInWord = BWTResidentSizeInWord(bwt->textLength) + WORD_BETWEEN_OCC / 2;
	bwtCodeLengthInFile = BWTFileSizeInWord(bwt->textLength);
	tempBwtCode = (unsigned int *)malloc(bwt->bwtSizeInWord * sizeof(unsigned int));
	printf("reference length: %u\n", bwt->textLength);	


	fread(tempBwtCode, sizeof(unsigned int), bwtCodeLengthInFile, bwtCodeFile);
	fclose(bwtCodeFile);

	DecodeBWT(bwt, tempBwtCode);
	BWTClearTrailingBwtCode(bwt);

	LoadLookupTable(bwt->lookupTable, lookupTableFileName, TABLE_SIZE);

#if LAYER_INDEX_2
	for(i=0; i<ALPHABET_SIZE+1; i++)
		fread(&tmp,sizeof(unsigned int), 1, occValueFile);
	fread(bwt->occValue, sizeof(unsigned int), bwt->occSizeInWord, occValueFile);
	fread(bwt->occMajorValue, sizeof(unsigned int), bwt->occMajorSizeInWord, occValueFile);
	fclose(occValueFile);
#endif

#if LAYER_INDEX_1
	uint32 temp;
	fread(&temp, sizeof(unsigned int), 1, gpuOccValueFile);
	for(i=0; i<ALPHABET_SIZE; i++){
		fread(&temp, sizeof(unsigned int), 1, gpuOccValueFile);
	}
	unsigned int gpuOccValueSize = (bwt->textLength + GPU_OCC_INTERVAL - 1) / GPU_OCC_INTERVAL * ALPHABET_SIZE;	
	bwt->GPUOccValue = (unsigned int *)malloc(gpuOccValueSize * sizeof(int));
	fread(bwt->GPUOccValue, sizeof(int), gpuOccValueSize, gpuOccValueFile);
	fclose(gpuOccValueFile);
	bwt->GPUOccValueInWord = gpuOccValueSize;

#endif

	revbwt->countValue = (unsigned int *)malloc((ALPHABET_SIZE + 1)*sizeof(int));
	fread(&revbwt->inverseSa0, sizeof(unsigned int), 1, revbwtCodeFile);
	fread(&revbwt->countValue[1], sizeof(unsigned int), ALPHABET_SIZE, revbwtCodeFile);
	Initialize(revbwt, revbwt->countValue[ALPHABET_SIZE]);
	revbwt->countValue[0] = 0;

	bwtCodeLengthInFile = BWTFileSizeInWord(revbwt->textLength);
	revbwt->bwtSizeInWord = BWTResidentSizeInWord(revbwt->textLength) + WORD_BETWEEN_OCC / 2;
	tempBwtCode = (unsigned int *)malloc(bwt->bwtSizeInWord * sizeof(unsigned int));
	fread(tempBwtCode, sizeof(unsigned int), bwtCodeLengthInFile, revbwtCodeFile);
	fclose(revbwtCodeFile);

	DecodeBWT(revbwt, tempBwtCode);
	BWTClearTrailingBwtCode(revbwt);
	LoadLookupTable(revbwt->lookupTable,revLookupTableFileName, TABLE_SIZE);

#if  LAYER_INDEX_1	
	fread(&temp, sizeof(unsigned int), 1, gpurevOccValueFile);
	for(i=0; i<ALPHABET_SIZE; i++)
		fread(&temp, sizeof(unsigned int), 1, gpurevOccValueFile);
	revbwt->GPUOccValue = (unsigned int *)malloc(gpuOccValueSize * sizeof(int));
	fread(revbwt->GPUOccValue, sizeof(int), gpuOccValueSize, gpurevOccValueFile);
	fclose(gpurevOccValueFile);
	revbwt->GPUOccValueInWord = gpuOccValueSize;
#endif

#if LAYER_INDEX_2
	for(i=0; i<ALPHABET_SIZE+1; i++)
		fread(&tmp,sizeof(unsigned int), 1,revoccValueFile);
	fread(revbwt->occValue, sizeof(unsigned int), revbwt->occSizeInWord, revoccValueFile);
	fread(revbwt->occMajorValue, sizeof(unsigned int), revbwt->occMajorSizeInWord, revoccValueFile);
	fclose(revoccValueFile);
#endif
}

void compressRead2Buffer(char* read, unsigned int* Queries, unsigned int readId, int queryLen)
{
	unsigned int i = readId;
	int j,k,pos = 0;
	int querySizeInWord = QUERY_SIZE_IN_WORD;

	for(j=0; j<querySizeInWord; j++) {
		for(k=0; k<CHAR_PER_WORD && pos<queryLen; k++) {
			Queries[i * querySizeInWord + j] <<= 2;
			Queries[i * querySizeInWord + j] |= (EncodeInt(read[pos]) & 0x3);
			pos++;
		}
		if(pos == queryLen) {
			for(; k<CHAR_PER_WORD; k++)
				Queries[i * querySizeInWord + j] <<= 2;
		}
	}
}

void addName2Buffer(char* readName, char* readNameBuffer, unsigned int* offset)
{
	int i = 0;
	while(readName[i] != '\0' && readName[i] != '\n'){
		readNameBuffer[(*offset)++] = readName[i++];
	}
	readNameBuffer[(*offset)++] = '\0';
}


bool isFASTQ(char * inputFileName) {
	bool FASTQ = false;

	if(strstr(inputFileName, "fastq") != NULL){
		FASTQ = true;
		printf("Input format: FASTQ.\n");
	} else {
		printf("Input format: FASTA.\n");
	}

	return FASTQ;
}

void parseInputFile(const int QUERY_NUM, FILE* fin, bool FASTQ, int* ReadLength, char* ReadNameBuffer, 
		unsigned int* ReadNum, unsigned int* Queries, 
		const bool reverse = false) {
	char c;
	unsigned int i=0, j=0, k=0, pos = 0;
	const int maxReadNameLength = 200;
	const int maxReadLength = 200;
	char ReadName[maxReadNameLength];
	char Read[maxReadLength];
	unsigned int ReadNameOffset = 0;
	char temp[200];


	fgets(ReadName, maxReadNameLength, fin);
	addName2Buffer(ReadName,ReadNameBuffer, &ReadNameOffset);
	fgets(Read, maxReadLength, fin);

	for(i=0; i<maxReadLength; i++) {
		if(Read[i] == '\n' || Read[i] == '\0')
			break;
	}
	(*ReadLength) = i;
	printf("%s query len: %d\n",Read, *ReadLength);

	int readLen = *ReadLength;
	for(j=0; j<(readLen); j++){
		c = Read[j];
		if(c == 'a' || c =='c' || c=='g' || c=='t'){
			Read[j] = toupper(c);
		}
		else if( c == 'A' || c =='C' || c=='G' || c=='T')
			Read[j] = c;
		else if(c == 'N' || c=='n')
			Read[j] = 'G';
	}
	unsigned int querySizeInWord = QUERY_SIZE_IN_WORD;
	i = 0;
	compressRead2Buffer(Read, Queries, i, readLen);
	i++;
	char* revQuery = (char*)malloc(sizeof(char)*(readLen));

	while(i < QUERY_NUM && !feof(fin))
	{
		fgets(ReadName,maxReadNameLength,fin);
		addName2Buffer(ReadName, ReadNameBuffer, &ReadNameOffset);
		//	printf("ReadNameBuffer: %s\n", ReadNameBuffer);
		//	printf("ReadNameOffset: %u\n", ReadNameOffset);

		fgets(Read, maxReadLength, fin);

		if(strlen(Read) != readLen + 1) {
			printf("%d different length ignore\n", strlen(Read));
			continue;
		}
		if(reverse == true) {
			reverse_complement(revQuery, Read, readLen);
			memcpy(Read, revQuery, sizeof(char)*readLen);
		}
		for(j=0; j<(readLen); j++){
			c = Read[j];
			if(c == 'a' || c =='c' || c=='g' || c=='t'){
				Read[j] = toupper(c);
			}
			else if( c == 'A' || c =='C' || c=='G' || c=='T')
				Read[j] = c;
			else if(c == 'N' || c=='n')
				Read[j] = 'G';
		}	

		compressRead2Buffer(Read,Queries, i, readLen);

		//pack Reads[i].query to Queries 
		if(FASTQ) {
			fgets(temp,200, fin);
			fgets(temp,200, fin);
		}
		i++;
	}
	(*ReadNum) = i;
	free(revQuery);				
}



//generate the reverse reads as well
void parseInputFileTwoStrand(const int QUERY_NUM, FILE* fin, bool FASTQ, int* ReadLength, char* ReadNameBuffer, 
		unsigned int* ReadNum, unsigned int* Queries, char* qualityBuffer = NULL) {
	char c;
	unsigned int i=0, j=0, k=0, pos = 0;
	const int maxReadNameLength = 200;
	const int maxReadLength = 200;
	char ReadName[maxReadNameLength];
	char Read[maxReadLength];
	char quality[maxReadLength];
	unsigned int ReadNameOffset = 0;
	char temp[200];

	fgets(ReadName, maxReadNameLength, fin);
	addName2Buffer(ReadName,ReadNameBuffer, &ReadNameOffset);
	fgets(Read, maxReadLength, fin);

	for(i=0; i<maxReadLength; i++) {
		if(Read[i] == '\n' || Read[i] == '\0')
			break;
	}
	(*ReadLength) = i;
	printf("%s query len: %d\n",Read, *ReadLength);
	int readLen = *ReadLength;
	char* revQuery = (char*)malloc(sizeof(char)*(readLen));

	if(FASTQ) {
		fgets(ReadName, maxReadNameLength, fin);
		fgets(qualityBuffer, maxReadLength, fin);
	}

	reverse_complement(revQuery, Read, readLen);

	for(j=0; j<(readLen); j++){
		c = Read[j];
		if(c == 'a' || c =='c' || c=='g' || c=='t'){
			Read[j] = toupper(c);
		}
		//else if( c == 'A' || c =='C' || c=='G' || c=='T')
		//	Read[j] = c;
		else if(c == 'N' || c=='n')
			Read[j] = 'G';
	}
	for(j=0; j<(readLen); j++){
		c = revQuery[j];
		if(c == 'a' || c =='c' || c=='g' || c=='t'){
			revQuery[j] = toupper(c);
		}
		//else if( c == 'A' || c =='C' || c=='G' || c=='T')
		//        Read[j] = c;
		else if(c == 'N' || c=='n')
			revQuery[j] = 'G';
	}
	unsigned int querySizeInWord = QUERY_SIZE_IN_WORD;
	i = 0;
	compressRead2Buffer(Read, Queries, i, readLen);
	compressRead2Buffer(revQuery, Queries, i + QUERY_NUM, readLen);


	i++;

	while(i < QUERY_NUM && !feof(fin))
	{
		fgets(ReadName,maxReadNameLength,fin);
		addName2Buffer(ReadName, ReadNameBuffer, &ReadNameOffset);

		fgets(Read, maxReadLength, fin);

		if(strlen(Read) != readLen + 1) {
			printf("%d different length ignore\n", strlen(Read));
			continue;
		}
		//if(reverse == true) {
		reverse_complement(revQuery, Read, readLen);
		//	memcpy(Read, revQuery, sizeof(char)*readLen);
		//}
		for(j=0; j<(readLen); j++){
			c = Read[j];
			if(c == 'a' || c =='c' || c=='g' || c=='t'){
				Read[j] = toupper(c);
			}
			//else if( c == 'A' || c =='C' || c=='G' || c=='T')
			//Read[j] = c;
			else if(c == 'N' || c=='n')
				Read[j] = 'G';
		}	


		for(j=0; j<(readLen); j++){
			c = revQuery[j];
			if(c == 'a' || c =='c' || c=='g' || c=='t'){
				revQuery[j] = toupper(c);
			}
			//else if( c == 'A' || c =='C' || c=='G' || c=='T')
			//	Read[j] = c;
			else if(c == 'N' || c=='n')
				revQuery[j] = 'G';
		}	

		compressRead2Buffer(Read, Queries, i, readLen);
		compressRead2Buffer(revQuery, Queries, i + QUERY_NUM, readLen);

		//pack Reads[i].query to Queries 
		
		if(FASTQ) {
			fgets(ReadName, maxReadNameLength, fin);
			fgets(qualityBuffer + i*readLen, maxReadLength, fin);
		}
/*		if(FASTQ) {
			fgets(temp,200, fin);
			fgets(temp,200, fin);
		}
*/
		i++;
	}
	(*ReadNum) = i;

	cout << "QUERY_NUM = " << QUERY_NUM << endl;
	cout << "*ReadNum = " << *ReadNum << endl; 

	memcpy(qualityBuffer +  (*ReadNum) * readLen, qualityBuffer, (*ReadNum) * readLen);
	assert(QUERY_NUM == *ReadNum);

	free(revQuery);				
}


void parseInputFile(const int QUERY_NUM, char * inputFileName, int* ReadLength, char* ReadNameBuffer, 
		unsigned int* ReadNum, unsigned int* Queries, 
		const bool reverse = false)
{
	char c;
	unsigned int i=0, j=0, k=0, pos = 0;
	const int maxReadNameLength = 200;
	const int maxReadLength = 200;
	char ReadName[maxReadNameLength];
	char Read[maxReadLength];
	FILE * fin = fopen(inputFileName, "r");
	unsigned int ReadNameOffset = 0;
	if(fin == NULL) {
		printf("can't open input file: %d\n", inputFileName);
	}
	char temp[200];
	bool FASTQ = false;

	if(strstr(inputFileName, "fastq") != NULL){
		FASTQ = true;
		printf("fastq\n");
	}

	fgets(ReadName, maxReadNameLength, fin);
	addName2Buffer(ReadName,ReadNameBuffer, &ReadNameOffset);
	fgets(Read, maxReadLength, fin);

	for(i=0; i<maxReadLength; i++) {
		if(Read[i] == '\n' || Read[i] == '\0')
			break;
	}
	(*ReadLength) = i;

	int readLen = *ReadLength;
	for(j=0; j<(readLen); j++){
		c = Read[j];
		if(c == 'a' || c =='c' || c=='g' || c=='t'){
			Read[j] = toupper(c);
		}
		else if( c == 'A' || c =='C' || c=='G' || c=='T')
			Read[j] = c;
		else if(c == 'N' || c=='n')
			Read[j] = 'G';
	}
	unsigned int querySizeInWord = QUERY_SIZE_IN_WORD;
	i = 0;
	compressRead2Buffer(Read, Queries, i, readLen);
	i++;
	char* revQuery = (char*)malloc(sizeof(char)*(readLen));

	while(i < QUERY_NUM && !feof(fin))
	{
		fgets(ReadName,maxReadNameLength,fin);
		addName2Buffer(ReadName, ReadNameBuffer, &ReadNameOffset);
		//	printf("ReadNameBuffer: %s\n", ReadNameBuffer);
		//	printf("ReadNameOffset: %u\n", ReadNameOffset);

		fgets(Read, maxReadLength, fin);

		if(strlen(Read) != readLen + 1) {
			printf("%d different length ignore\n", strlen(Read));
			continue;
		}
		if(reverse == true) {
			reverse_complement(revQuery, Read, readLen);
			memcpy(Read, revQuery, sizeof(char)*readLen);
		}
		for(j=0; j<(readLen); j++){
			c = Read[j];
			if(c == 'a' || c =='c' || c=='g' || c=='t'){
				Read[j] = toupper(c);
			}
			else if( c == 'A' || c =='C' || c=='G' || c=='T')
				Read[j] = c;
			else if(c == 'N' || c=='n')
				Read[j] = 'G';
		}	

		compressRead2Buffer(Read,Queries, i, readLen);

		//pack Reads[i].query to Queries 
		if(FASTQ) {
			fgets(temp,200, fin);
			fgets(temp,200, fin);
		}
		i++;
	}
	(*ReadNum) = i;
	fclose(fin);	
	free(revQuery);
}


uint64 getNumAlign(uint3* sa, const uint32 numSA) {
	uint64 numAlign = 0;
	for(uint64 i = 0; i < numSA; i++) {
		numAlign += (sa[i].z - sa[i].y + 1);
	}

	return numAlign;
}

void sortQuery(uint3* cpuSA, unsigned int* Queries, const int numQuery, 
		BWT* bwt, BWT* revbwt) {
	preprocess(cpuSA, Queries, numQuery, bwt, revbwt);

	FILE* file = NULL;
	safe_fopen(&file, "query.bin", "wb");
	safe_fwrite(Queries, sizeof(unsigned int), QUERY_SIZE_IN_WORD*numQuery, file);
	safe_fclose(file);
}


void parseBinInputFile(const char* fileName, unsigned int* Queries, const int numQuery) {
	FILE* file = NULL;
	safe_fopen(&file, fileName, "rb");
	safe_fread(Queries, sizeof(unsigned int), QUERY_SIZE_IN_WORD*numQuery, file);
	safe_fclose(file);
}


void calculateNumMatch(long* numMatch, const uint3* sa, const uint32 numSA) {
	for(int i = 0; i < numSA; i++) {
		numMatch[sa[i].x] += (sa[i].z - sa[i].y + 1);
	}
}
/*
uint32 getRefPos(uint3* SA_result, uint32 numSA, uint32* saIndex, uint32* queryId,uint32* offset, uint32* pos)
{
	uint32 i,j, k;
	k = 0;
	printf("numSA: %u\n", numSA);
	for(i=0; i<numSA; i++){
		queryId[i] = SA_result[i].x;
		offset[i] = k;
		for(j=SA_result[i].y; j<=SA_result[i].z; j++) {
			if(j > 2858012806) {
				printf("%u\n", j);
				continue;
			}
			pos[k] = saIndex[j];
			k++;
		}
	}
	printf("num alignment: %u\n", k);
	return k;
}
*/

struct Index {
	//for CPU
	BWT* bwt;
	BWT* revbwt;
	uint32* saIndex;

	//for GPU
	BWT* d_bwt;
	BWT* d_revbwt;
	st_MemoryBuf* GPUMemBuf;
};

/*
struct Buffer {
	//for alignment result output
	uint64 saBufSize;
	uint3* cpuSA;

	//for getRefPos
	//these three pointers shoule be allocated using allocate after the alignment process
	uint32* queryId;
	uint32* offset;
	uint32* pos;
	//these two pointers are allocated
	uint32* c_queryId;
	uint32* c_queryId_len;

	void allocate(const int numSA, const int nresult) {
		queryId =  (uint32*)malloc(sizeof(int)*numSA);
		offset =  (uint32*)malloc(sizeof(int)*numSA);
		pos =  (uint32*)malloc(sizeof(int)*nresult);

		printf("Buffer::allocate done: %lu bytes.\n", (uint64)sizeof(int)*(numSA*2 + (uint64)nresult));
	}

	void release() {
		free(queryId);
		free(offset);
		free(pos);
	}
};


struct Read{
	char* readNameBuffer;
	unsigned int* readBuffer;
	uint32 numQuery;

	int textLength;
	int readNum;
	int maxReadNameLength;
	int readLength; 

	void writeMeta(FILE* resultFile) {
		fwrite(&textLength, sizeof(int), 1, resultFile);
		fwrite(&readNum, sizeof(int), 1, resultFile);
		fwrite(&maxReadNameLength, sizeof(int), 1, resultFile);
		fwrite(&readLength, sizeof(int), 1, resultFile);
	}
};

*/



void runSinglePassNoCompression(FILE* resultFile, Read* read, Index* index, Buffer* buf, AlignResult* alignRes, bool compressOutput = true, bool limited = false, char* ref = NULL) {

	unsigned int timer = 0;
	uint32 nresult = 0;
	/*allocate the CPU result buffer*/
	//buf->saBufSize = 8 *((uint64)(1024*1024*1024));
	//buf->cpuSA = (uint3*)malloc(buf->saBufSize);

	/*alignment*/
	uint32 numSA = 0;
	const int numCPUThread = 16;
	index->GPUMemBuf = new st_MemoryBuf(GPU_BUFSIZE, GMEM);
	startTimer(&timer);
	//GPU-CPU coprocessing for all alignments or for limit

	/////////////////
	//nresult = coproc2BWTSearch_simple(read->readBuffer, read->numQuery, 
	//			index->bwt, index->revbwt, index->d_bwt, index->d_revbwt, index->GPUMemBuf,
	//			numCPUThread, buf->saBufSize, buf->cpuSA, &numSA);
	//const int nresult = coproc2BWTSearchLimit_simple(read->readBuffer, read->numQuery, 
	//                        index->bwt, index->revbwt, index->d_bwt, index->d_revbwt, index->GPUMemBuf,
	//                        numCPUThread, buf->saBufSize, buf->cpuSA, &numSA);
	/////////////////

	/*GPU for all begins*/
	nresult = GPU2BWTSearch_shareBuf(read->readBuffer, read->numQuery, index->d_bwt, index->d_revbwt,
                			index->GPUMemBuf, buf->cpuSA, &numSA);
	/*GPU for all ends*/
	/*GPU for limit	begins */
	////////////////////
	/*
	if(limited) {
		printf("numQuery = %d\n", read->numQuery);
		long* d_numAlign = NULL;
		GPUMALLOC((void**)&d_numAlign, sizeof(long)*read->numQuery);
		cutilSafeCall(cudaMemset(d_numAlign, 0, sizeof(long)*read->numQuery));
		GPU2BWTSearchLimit(d_numAlign, read->readBuffer, read->numQuery, 
			index->d_bwt, index->d_revbwt, index->GPUMemBuf, 
			buf->cpuSA, &numSA);
		long* gpu_numAlign = (long*)malloc(sizeof(long)*read->numQuery );
		FROMGPU(gpu_numAlign, d_numAlign, sizeof(long)*read->numQuery);
		nresult = reportNumMatch(gpu_numAlign, read->numQuery);

		printf("\n\n");
		free(gpu_numAlign);
		GPUFREE(d_numAlign);
	}*/
	///////////////////
	/*GPU for limit ends*/
	cout << "numAlign = " << nresult << endl;
	cudaThreadSynchronize();
	alignSec += endTimer(&timer);
	index->GPUMemBuf->close();


	/*allocate temp buffer*/
	buf->allocate(numSA, nresult);

	/*compression*/
	timer = 0;
	startTimer(&timer);

	uint32 outputSize = 0;
	sortSA(buf->cpuSA, numSA);
	printf("===> sortSA done.\n");
	getRefPos(buf->cpuSA, numSA, index->saIndex, nresult, buf->queryId, buf->offset, buf->pos);
	printf("===> getRefPos done.\n");


////////////////////////////////////////
//debug only
//	sort(buf->pos, buf->pos + nresult);	
	/*for(int i = 0; i < 100; i++) {
		cout << buf->pos[i] << endl;
	}*/

/////////////////////////////////////


	/*allocate memory for alignRes*/
	alignRes->alloc(nresult, read->readLength);
	cout << "alignRes->_readLength = " << alignRes->_readLength << endl;
	cout << "alignRes->_numAlign = " << alignRes->_numAlign << endl;

	//	fullSizeOutput(read->readNameBuffer, read->readBuffer, buf);
	//findInfo(read, buf, nresult, numSA, alignRes, ref);
/////////////////////////////////////////////
//for temp debug, you should find the complete information later
	memcpy(alignRes->_pos, buf->pos, sizeof(unsigned int)*nresult);
//////////////////////////////////////////	


	/*free temp buffer*/
	buf->release();
	//free(buf->cpuSA);
	printf("free memory done.\n");
}

void runSinglePass(FILE* resultFile, Read* read, Index* index, Buffer* buf, bool compressOutput = true, bool limited = false, char* ref = NULL) {

	unsigned int timer = 0;
	uint32 nresult = 0;
	/*allocate the CPU result buffer*/
	//buf->saBufSize = 8 *((uint64)(1024*1024*1024));
	//buf->cpuSA = (uint3*)malloc(buf->saBufSize);




	/*alignment*/
	uint32 numSA = 0;
	const int numCPUThread = 16;
	index->GPUMemBuf = new st_MemoryBuf(GPU_BUFSIZE, GMEM);
	startTimer(&timer);
	//GPU-CPU coprocessing for all alignments or for limit

	/////////////////
	//nresult = coproc2BWTSearch_simple(read->readBuffer, read->numQuery, 
	//			index->bwt, index->revbwt, index->d_bwt, index->d_revbwt, index->GPUMemBuf,
	//			numCPUThread, buf->saBufSize, buf->cpuSA, &numSA);
	//const int nresult = coproc2BWTSearchLimit_simple(read->readBuffer, read->numQuery, 
	//                        index->bwt, index->revbwt, index->d_bwt, index->d_revbwt, index->GPUMemBuf,
	//                        numCPUThread, buf->saBufSize, buf->cpuSA, &numSA);
	/////////////////



	/*GPU for all begins*/
	nresult = GPU2BWTSearch_shareBuf(read->readBuffer, read->numQuery, index->d_bwt, index->d_revbwt,
                			index->GPUMemBuf, buf->cpuSA, &numSA);
	/*GPU for all ends*/


	
	/*GPU for limit	begins */
	////////////////////
	/*
	if(limited) {
		printf("numQuery = %d\n", read->numQuery);
		long* d_numAlign = NULL;
		GPUMALLOC((void**)&d_numAlign, sizeof(long)*read->numQuery);
		cutilSafeCall(cudaMemset(d_numAlign, 0, sizeof(long)*read->numQuery));
		GPU2BWTSearchLimit(d_numAlign, read->readBuffer, read->numQuery, 
			index->d_bwt, index->d_revbwt, index->GPUMemBuf, 
			buf->cpuSA, &numSA);
		long* gpu_numAlign = (long*)malloc(sizeof(long)*read->numQuery );
		FROMGPU(gpu_numAlign, d_numAlign, sizeof(long)*read->numQuery);
		nresult = reportNumMatch(gpu_numAlign, read->numQuery);

		printf("\n\n");
		free(gpu_numAlign);
		GPUFREE(d_numAlign);
	}*/
	///////////////////
	/*GPU for limit ends*/

	cudaThreadSynchronize();
	alignSec += endTimer(&timer);
	index->GPUMemBuf->close();


	/*allocate temp buffer*/
	buf->allocate(numSA, nresult);

	/*compression*/
	timer = 0;
	startTimer(&timer);
	uint32 compressedSize = 0;
	uint32 compressedTitleSize = 0;
	uint32 outputSize = 0;
	sortSA(buf->cpuSA, numSA);
	printf("===> sortSA done.\n");
	getRefPos(buf->cpuSA, numSA, index->saIndex, nresult, buf->queryId, buf->offset, buf->pos);
	printf("===> getRefPos done.\n");

	//	fullSizeOutput(read->readNameBuffer, read->readBuffer, buf);
	if(limited) {
		read->writeMeta(resultFile);
		compressTitle(read->readNameBuffer, read->readNum, resultFile, &compressedTitleSize);
		printf("===> compressTitle done.\n");
		compressQueryId(buf->queryId, numSA, read->numQuery, buf->c_queryId, buf->c_queryId_len, &compressedSize);
		printf("===> compressQueryId done.\n");
		outputIdPos(buf->c_queryId, buf->c_queryId_len, compressedSize, buf->offset, 
				numSA, buf->pos, nresult, resultFile, &outputSize);
		printf("===> outputIdPos done.\n");
		compressOutputLimited(read, buf, nresult, numSA, resultFile, ref);
		printf("===> compress&output mis info done.\n");
		cudaThreadSynchronize();	
	}
	else if(!compressOutput) {
		fullSizeOutput(read,buf,nresult,numSA, resultFile, ref);
	}
	else {
		read->writeMeta(resultFile);
		compressTitle(read->readNameBuffer, read->readNum, resultFile, &compressedTitleSize);
		printf("===> compressTitle done.\n");
		fwrite(read->readBuffer, sizeof(int), QUERY_SIZE_IN_WORD*read->readNum, resultFile);
		compressQueryId(buf->queryId, numSA, read->numQuery, buf->c_queryId, buf->c_queryId_len, &compressedSize);
		printf("===> compressQueryId done.\n");
		cudaThreadSynchronize();
		comprSec += endTimer(&timer);		

		/*dump to the disk*/
		timer = 0;
		startTimer(&timer);
		outputIdPos(buf->c_queryId, buf->c_queryId_len, compressedSize, buf->offset, 
				numSA, buf->pos, nresult, resultFile, &outputSize);
		printf("===> outputIdPos done.\n");
		cudaThreadSynchronize();
		outputSec += endTimer(&timer);
	}
	/*free temp buffer*/
	buf->release();
	//free(buf->cpuSA);
	printf("free memory done.\n");
}
void loadRef(char* ref, uint32 textLength, char* indexFileName)
{
	char refFileName[200];
	int refFileNameLength = strstr(indexFileName, ".index") - indexFileName;
	strncpy(refFileName, indexFileName, refFileNameLength);
	refFileName[refFileNameLength] = '\0';
	FILE * refFile = fopen(refFileName, "r");
	char c;
	uint32 i;
	c = fgetc(refFile);
	if(c != '>') {
		ref[0] = c;
		i = 1;
	}
	else {
		while(c != '\n'){
			c = fgetc(refFile);
		}
		i = 0;
	}
	
	while(i < textLength && !feof(refFile)) {
		ref[i] = fgetc(refFile);
		
		if(ref[i] != '\n' && ref[i] != 'n' && ref[i] != 'N') {
			ref[i] = toupper(ref[i]);
		}
/*	
		int temp = 0;
		while((ref[i] == 'N' || ref[i] == 'n') && temp < 10) {
			temp++;
			ref[i] = 'G';
			i++;	
			ref[i] = fgetc(refFile);		
		}	
	
		if(temp == 10) {
			i = i-10;
			while(ref[i] == 'N' || ref[i] == 'n') {
				ref[i] = fgetc(refFile);
			}
		}

*/		if(ref[i] != '\n')
			i++;
	}
	printf("load ref done length = %u\n", i);	
	fclose(refFile);
	
}

/*
TODO:
1. parse the input read file one time.
2. due to assert(QUERY_NUM == *ReadNum);, currntly only support each pass reads the same #reads
 */
int main(int argc, char* argv[]) {

	unsigned int timer = 0;
	unsigned int totalTimer = 0;
	bool compress = true, limited = false;
	BWT *bwt, *revbwt;
	char *ref = NULL;
	startTimer(&timer);
	//load the bwt index
	char* indexFileName = argv[1];
	char* queryFileName = argv[2];
	char* resultFileName = argv[3];

	for(int i=4; i<argc; i++) {
		if(strcmp(argv[i], "-nc") == 0)
			compress = false;
		if(strcmp(argv[i], "-l") == 0)
			limited = true;
	}
	
	bwt = (BWT *)malloc(sizeof(BWT));
	bwt->lookupTable = (LookupTable *)malloc(sizeof(LookupTable));
	revbwt = (BWT *)malloc(sizeof(BWT));
	revbwt->lookupTable = (LookupTable *)malloc(sizeof(LookupTable));
	LoadIndex(bwt, revbwt, indexFileName);
	printf("Load index done. \n");

	//load the SA index
	uint32 textLength = bwt->textLength;
	if(compress == false || limited == true) {
		ref = (char*)malloc(textLength);
		loadRef(ref, textLength, indexFileName);
	}

	uint32* saIndex = NULL;
	char saFileName[200] = "";
	const char* saFileSuffix = ".sa";
	strcat(saFileName, indexFileName);
	strcat(saFileName, saFileSuffix);
	FILE* saFile = NULL;
	safe_fopen(&saFile, saFileName, "rb");
	uint32 saInterval;
	for(int i=0; i<6; i++) {
		fread(&saInterval, sizeof(int), 1, saFile);
	}
	
	uint32 saIndexSize = (textLength + saInterval) / saInterval;
	saIndex = (uint32*)malloc(saIndexSize * sizeof(int));
	safe_fread(saIndex, sizeof(int), saIndexSize, saFile);
	safe_fclose(saFile);
	printf("Load SA done.\n");

	cudaThreadSynchronize();
	endTimer(&timer, "All index loading");

	printf("\nNow index data ready!\n");
	printf("-------------------------------------------------------\n\n\n");


	cudaThreadSynchronize();
	startTimer(&totalTimer);

	uint32 numSA1 = 0, compressedTitleSize = 0;
	uint64 numAlign1, numAlign2;
	uint64 totalOutputSize = 0;
	int readLength;



	const int numTotalQuery = 1000000;
	const int numQueryPerPass = 1000000;
	assert(numTotalQuery%numQueryPerPass == 0);


	//init memory
	int i,j;
	int matchCounter = 0;
	unsigned int matchRead = 0;
	unsigned int ReadNum = 0;
	const int maxReadNameLength = 200;
	const int maxReadLength = 100;
	const int maxReadLengthInWord = maxReadLength / 16;
	char* readNameBuffer = (char*)malloc(numQueryPerPass * maxReadNameLength);
	uint32* readBuffer = (uint32*)malloc(numQueryPerPass * 2 * sizeof(int) * maxReadLengthInWord);
	char* qualityBuffer = (char*)malloc(numQueryPerPass * 2 * maxReadLength);

	//initialize the cpu result buffer
	//uint64 saBufSize = 8 *((uint64)(1024*1024*1024));
	//uint3* cpuSA = (uint3*)malloc(saBufSize);
	Read* read = (Read*)malloc(sizeof(Read));
	Index* index = (Index*)malloc(sizeof(Index));
	Buffer* buffer = (Buffer*)malloc(sizeof(Buffer));
	AlignResult * alignRes = new AlignResult;
	
	
	buffer->c_queryId = (uint32*)malloc(sizeof(uint32)*numQueryPerPass*2);
	buffer->c_queryId_len = (uint32*)malloc(sizeof(uint32)*numQueryPerPass*2);



	//open the result file
	FILE* resultFile = NULL;
	safe_fopen(&resultFile, resultFileName, "wb");


	//init the CPU lookup tables
	CPU2BWTInit(bwt, revbwt);

	//init the GPU memory
	const uint64 GPUBufSize = (2.8*((uint64)(1024*1024*1024)));
	GPU2BWTInitIndex(&index->d_bwt, &index->d_revbwt, bwt, revbwt, &index->GPUMemBuf, GPUBufSize);

	//open the input file	
	FILE * fin = NULL;
	safe_fopen(&fin, queryFileName, "r");
	int querySizeInWord = QUERY_SIZE_IN_WORD;	
	bool fastq = isFASTQ(queryFileName);


	/*allocate the CPU result buffer*/
	buffer->saBufSize = 4 *((uint64)(1024*1024*1024));
	buffer->cpuSA = (uint3*)malloc(buffer->saBufSize);	

	/*prepare the output container*/
	const unsigned int refLength = 3000000000;
	const unsigned int winSize = 32768;
	const unsigned int bucketSize = 100;
	string outFileName = "./result.bin";
	AlignResultWrite* alignResultWrite = new AlignResultWrite(QUERY_LENGTH, 
					refLength, winSize, bucketSize, outFileName);

	const int numPass = numTotalQuery/numQueryPerPass;
	assert(numTotalQuery%numQueryPerPass == 0);
	for(int passId = 0; passId < numPass; passId++) {

		printf("pass %d ........................\n", passId);

		timer = 0;
		startTimer(&timer);
		parseInputFileTwoStrand(numQueryPerPass, fin, fastq, &readLength, 
				readNameBuffer, &ReadNum, readBuffer, qualityBuffer);
		assert(ReadNum == numQueryPerPass);
		inputSec += endTimer(&timer);

		const uint32 numQuery = ReadNum*2;
		read->readNameBuffer = readNameBuffer;
		read->readBuffer = readBuffer;
		read->numQuery = numQuery;
		read->textLength = textLength;
		read->readNum = ReadNum;
		read->maxReadNameLength = maxReadNameLength;
		read->readLength = readLength;
		index->bwt = bwt;
		index->revbwt = revbwt;
		index->saIndex = saIndex;
		//buffer->saBufSize = saBufSize;
		//buffer->cpuSA = cpuSA;
		assert(read->readLength == QUERY_LENGTH);


		//runSinglePass(resultFile, read, index, buffer, compress, limited, ref);
		runSinglePassNoCompression(resultFile, read, index, buffer, alignRes, compress, limited, ref);

		//insert to the AlignResultWrite and release the memory
		alignRes->_qscore = qualityBuffer;
		alignResultWrite->insert(alignRes);
		alignRes->clear();
	}

	alignResultWrite->close();

	free(buffer->c_queryId);
	free(buffer->c_queryId_len);
	free(buffer->cpuSA);
	safe_fclose(fin);
	delete alignResultWrite;

	printf("\n\n ------------------------------------- \n");
	printf("All work done!\n");	
	printf("inputSec = %.3f sec\n", inputSec);
	printf("alginSec = %.3f sec\n", alignSec);
	printf("comprSec = %.3f sec\n", comprSec);
	printf("outputSec = %.3f sec\n", outputSec);

	cudaThreadSynchronize();
	endTimer(&totalTimer, "Overall processing time (exclude the index loading)");
}

